#define CV_NO_BACKWARD_COMPATIBILITY

#include <cxcore.h>
#include <cv.h>
#include <highgui.h>
#include <vector>
#include <iostream>

// Rearrange the quadrants of Fourier image so that the origin is at
// the image center
// src & dst arrays of equal size & type
void shiftDFT(cv::Mat & src_arr, cv::Mat & dst_arr )
{
    cv::Size size = src_arr.size();
    cv::Size dst_size = dst_arr.size();

    if(dst_size.width != size.width || dst_size.height != size.height)
    {
        //cvError( CV_StsUnmatchedSizes, "cvShiftDFT", "Source and Destination arrays must have equal sizes", __FILE__, __LINE__ );
        return;
    }

    int cx = size.width/2;
    int cy = size.height/2; // image center

    cv::Mat q1(src_arr, cv::Rect(0, 0, cx, cy));
    cv::Mat q2(src_arr, cv::Rect(cx, 0, cx, cy));
    cv::Mat q3(src_arr, cv::Rect(cx, cy, cx, cy));
    cv::Mat q4(src_arr, cv::Rect(0, cy, cx, cy));
    cv::Mat d1(dst_arr, cv::Rect(0, 0, cx, cy));
    cv::Mat d2(dst_arr, cv::Rect(cx, 0, cx, cy));
    cv::Mat d3(dst_arr, cv::Rect(cx, cy, cx, cy));
    cv::Mat d4(dst_arr, cv::Rect(0, cy, cx, cy));

    if(src_arr.data!=dst_arr.data)
    {
        if(q1.type() != d1.type())
        {
            //cvError( CV_StsUnmatchedFormats, "cvShiftDFT", "Source and Destination arrays must have the same format", __FILE__, __LINE__ );
            return;
        }
        q1.copyTo(d3);
        q2.copyTo(d4);
        q3.copyTo(d1);
        q4.copyTo(d2);
    }
    else
    {
        cv::Mat tmp(size, src_arr.type());
        q1.copyTo(tmp);
        q3.copyTo(q1);
        tmp.copyTo(q3);
        q2.copyTo(tmp);
        q4.copyTo(q2);
        tmp.copyTo(q4);
    }
}

int main(int argc, char ** argv)
{
    const char* filename = argc >=2 ? argv[1] : "lena.jpg";

    cv::Mat im = cv::imread( filename, CV_LOAD_IMAGE_GRAYSCALE);
    if( !im.data)
        return -1;

    cv::Mat realInput(im.size(), CV_64F);
    im.convertTo(realInput, CV_64F, 1.0, 0.0);
    cv::Mat imaginaryInput = cv::Mat::zeros(im.size(), CV_64F);
    cv::Mat complexInput(im.size(), CV_64FC2);

    cv::vector<cv::Mat> v;  // not std::vector
    v.push_back(realInput);
    v.push_back(imaginaryInput);
    cv::merge(v, complexInput);

    cv::Mat im2(im.size(), CV_8U);
    realInput.convertTo(im2, CV_8U, 1.0, 0.0);
    ///cv::imshow("win1", im2);

    ////
    std::cout<<"----------------------------------\n";
    std::cout<<"im type: "<<im.type()<<std::endl;
    std::cout<<CV_64F<<" "<<CV_64FC1<<" "<<CV_8U<<" "<<CV_8UC1<<std::endl;
    std::cout<<static_cast<int>(im.at<unsigned char>(200, 200))<<" "
        <<realInput.at<double>(200, 200)<<std::endl;

    int dft_M = cv::getOptimalDFTSize( im.rows - 1 );
    int dft_N = cv::getOptimalDFTSize( im.cols - 1 );

    ////
    std::cout<<"----------------------------------\n";
    std::cout<<"dft_M, dft_N "<<dft_M<<" "<<dft_N<<std::endl; 

    cv::Mat dft_A(dft_M, dft_N, CV_64FC2);

    // copy A to dft_A and pad dft_A with zeros
    cv::Mat tmp(dft_A, cv::Rect(0, 0, im.rows, im.cols));
    complexInput.copyTo(tmp);
    if( dft_A.cols > im.cols )
    {
        cv::Mat tmp(dft_A, cv::Rect(im.cols, 0, dft_A.cols-im.cols, im.rows));
        tmp = cv::Scalar(0, 0);
    }

    // no need to pad bottom part of dft_A with zeros because of
    // use nonzero_rows parameter in cvDFT() call below

    cv::dft(dft_A, dft_A, CV_DXT_FORWARD, complexInput.rows);

    cv::imshow("win", im);

    cv::vector<cv::Mat> res;
    // Split Fourier in real and imaginary parts
    cv::split(dft_A, res);
    cv::Mat image_Re(res[0]);
    cv::Mat image_Im(res[1]);

    ///cv::imshow("magnitude0", image_Re);
    ///cv::imshow("magnitude1", image_Im);

    std::cout<<image_Re.at<double>(200, 200)<<std::endl;
    // Compute the magnitude of the spectrum Mag = sqrt(Re^2 + Im^2)
    cv::pow(image_Re, 2.0, image_Re);

    std::cout<<image_Re.at<double>(200, 200)<<std::endl;

    cv::pow(image_Im, 2.0, image_Im);
    cv::add(image_Re, image_Im, image_Re);
    cv::pow(image_Re, 0.5, image_Re);

    // Compute log(1 + Mag)
    cv::add(image_Re, 1.0, image_Re); // 1 + Mag
    cv::log(image_Re, image_Re); // log(1 + Mag)


    // Rearrange the quadrants of Fourier image so that the origin is at
    // the image center
    shiftDFT(image_Re, image_Re);

    double m, M;
    cv::minMaxLoc(image_Re, &m, &M);
    std::cout<<image_Re.at<double>(200, 200)<<std::endl;
    //cv::convertScaleAbs(image_Re, image_Re, 1.0/(M-m), 1.0*(-m)/(M-m));
    //image_Re.convertTo(image_Re, CV_64F, 1.0/(M-m), 1.0*(-m)/(M-m));
    //image_Re=(image_Re - cv::Scalar(m))*(1.0/(M-m));// cannot compile
#if 0
    for (int i=0; i<512; ++i)
        for (int j=0; j<512; ++j)
            image_Re.at<double>(i,j) = (image_Re.at<double>(i,j) -m)/(M-m);
#else
    image_Re.convertTo(image_Re, CV_64F, 1.0/(M-m), 1.0*(-m)/(M-m));
#endif
    std::cout<<image_Re.at<double>(200, 200)<<std::endl;
    std::cout<<"m, M "<<m<<" "<<M<<std::endl; 

    //image_Re.convertTo(im2, CV_8U, 1.0, 0.0); //not work
    //cv::imshow("im2", im2);

    cv::imshow("magnitude", image_Re);

    cv::waitKey(-1);
    return 0;
}
